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1  Introduction 


Breast  cancer  is  the  most  frequently  diagnosed  malignancy  among  women  in  the  United 
States  [1].  In  1996,  the  American  Cancer  Society  estimated  that  184,300  women  would  be 
newly  diagnosed  with  breast  cancer  and  that  44,300  would  die  from  the  disease  [1].  Breast 
cancer  accounts  for  31%  of  all  cancers  detected  and  17%  of  all  cancer  deaths,  and  ranks  as 
the  second  leading  cause  of  death  from  cancer  among  women  in  the  United  States  [1].  Five 
year  survival  rates  are  generally  very  high  (93%)  for  breast  cancer  staged  as  being 
localized,  falling  to  72%  for  regional  disease  and  only  18%  for  distant  disease  [2].  The  early 
detection  of  breast  cancer  is  clearly  a  key  ingredient  of  any  strategy  designed  to  reduce 
breast  cancer  mortality. 

Mammography’s  role  is  the  early  detection  of  breast  cancer.  Although  more  accurate  than 
any  other  modality,  existing  techniques  for  mammography  only  find  80  to  90%  of  the  breast 
cancers.  Moreover,  in  7  to  10%  of  cases,  the  cancer  will  not  be  visible  on  the  mammogram. 
It  has  been  suggested  that  mammograms  as  normally  viewed,  display  only  about  3%  of  the 
total  information  detected.  Perception  is  a  problem  particularly  for  patients  with  dense 
fibroglandular  patterns.  The  importance  of  diagnosis  of  breast  cancer  at  an  early  stage  is 
critical  to  patient  survival.  The  general  inability  to  detect  small  tumors  and  other  salient 
features  within  mammograms  motivates  our  investigation. 

The  goal  of  this  project  is  to  develop  a  diagnostic  tool  for  radiologists  that  will  refine  the 
perception  of  mammographic  features  (including  lesions,  masses  and  calcifications)  and 
improve  the  accuracy  of  diagnosis.  Our  research  efforts  are  geared  towards  improving  the 
local  mammographic  viewing  environment  by  selectively  processing  mammograms  for 
presence  of  different  features,  and  towards  providing  a  better  global  mammographic 
viewing  environment  by  fusing  together  locally  processed  sections  of  images.  By  improving 
the  visualization  of  breast  pathology,  we  can  increase  the  chances  of  early  detection  of 
breast  cancers  (improve  quality)  while  requiring  less  time  to  evaluate  mammograms  for 
most  patients  (lower  costs). 

A  major  reason  for  poor  visualization  of  small  malignant  masses  is  the  subtle  difference  in 
x-ray  attenuation  between  normal  glandular  tissues  and  malignant  disease  [3] .  This  fact 
makes  the  detection  of  small  malignancies  problematical,  especially  in  younger  women  who 
have  denser  breast  tissue.  Although  calcifications  have  high  inherent  attenuation 
properties,  their  small  size  also  results  in  a  low  subject  contrast  [4].  As  a  result,  the 
visibility  of  small  tumors,  and  any  associated  microcalcifications,  will  always  be  a  problem 
in  mammography  as  it  is  currently  performed  using  analog  film. 

We  are  investigating  a  methodology  for  accomplishing  mammographic  feature  analysis 
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through  multiscale  representations.  Wide  variety  of  feature  sizes  and  shapes  in 
mammograms  makes  single-scale  processing  methods  ineffective.  In  [5],  it  was  shown  that, 
in  the  context  of  mammography,  multiscale  image  processing  algorithms  can  outperform 
traditional  contrast  enhancement  methods  such  as  histogram  equalization  and  unsharp 
masking.  As  reported  in  [6],  an  improvement  in  feature  visualization  was  noted  for 
mammograms  processed  using  multiscale  wavelet  processing  techniques.  Furthermore,  in 
[7],  it  was  demonstrated  that  unsharp  masking  with  a  Gaussian  lowpass  filter  can  be 
formulated  as  a  special  case  of  contrast  enhancement  via  a  discrete  dyadic  wavelet 
transform. 

Here,  we  present  a  scheme  for  local  enhancement  and  fusion  of  clinically  significant 
features.  We  devised  a  wavelet  transform  that  is  flexible  enough  for  incorporation  of  a 
variety  of  enhancement  methods  and  used  the  derived  wavelet  framework  for  enhancement 
of  microcalcifications,  circumscribed  masses,  and  stellate  lesions. 

In  the  following  sections,  we  briefly  overview  the  contents  of  the  report,  list  publications, 
and  explain  notation  that  we  use. 

1.1  Overview  of  Contents 

Mammographic  image  enhancement  methods  are  typically  aimed  at  either  improvement  of 
the  overall  visibility  of  features  or  enhancement  of  a  specific  sign  of  malignancy.  In  this 
report,  we  present  a  synthesis  of  the  two  paradigms  by  means  of  image  fusion.  Section  2.1 
introduces  the  idea  and  provides  an  overview  of  the  algorithm. 

Enhancement  of  different  mammographic  features  is  achieved  within  a  single  wavelet 
transform  framework  which  must  enable  inclusion  of  a  variety  of  methods  (the  derivations 
of  such  a  versatile  transform  can  be  viewed  as  an  extension  of  Task  3  of  our  Statement  of 
Work).  The  transform  is  described  in  Section  2.2.  First,  Section  2.2.1  builds  upon  the 
foundation  from  our  previous  report  to  derive  the  transform,  and  then,  Section  2.2.2  deals 
with  fast  implementations  issues. 

Different  methods  are  used  to  enhance  distinct  mammographic  features.  Section  2.3 
reports  on  three  enhancement  strategies.  Enhancement  of  microcalcifications  is  based  upon 
second  derivatives  of  a  Gaussian  and  given  in  Section  2.3.1.  Circumscribed  masses  are 
enhanced  through  inputting  the  result  from  Laplacian  of  Gaussian  approximations  to 
enhancement  functions  as  shown  in  Section  2.3.2.  Stellate  lesions  are  the  topic  of  Section 
2.3.3:  directional  derivatives  are  employed  to  determine  local  orientations  needed  for 
enhancement. 

Lastly,  Section  2.4  mentions  the  framework  for  fusion  of  enhanced  features. 
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1.2  Publications 


Below,  we  provide  the  list  of  publications  accomplished  during  1998. 

[1]  I.  Koren,  A.  Laine,  and  F.  Taylor,  “Enhancement  via  fusion  of  mammographic 
features,”  in  Proc.  IEEE  Int.  Conf.  Image  Process.,  Chicago,  II,  Oct.  1998. 

[2]  I.  Koren,  A.  Laine,  S.  Smith,  E.  Nickoloff,  and  F.  Taylor,  “Enhancement  of 
mammograms  via  fusion  of  enhanced  features,”  in  1st  Int.  Workshop  on  Computer-Aided 
Diagnosis,  Chicago,  II,  Sep.  1998. 

[3]  I.  Koren  and  A.  Laine,  “A  discrete  dyadic  wavelet  transform  for  multidimensional 
feature  analysis,”  in  M.  Akay  (Editor),  Time- Frequency  and  Wavelets  in  Biomedical  Signal 
Engineering,  New  York,  NY:  IEEE  Press,  1998,  pp.  425-449. 

1.3  Notation 

We  use  symbols  N,  Z,  and  R  for  the  sets  of  naturals,  integers,  and  reals,  respectively. 
L2{R )  and  L2(R2)  denote  the  Hilbert  spaces  of  measurable,  square-integrable  functions 
f(x )  and  f(x,y),  respectively. 

The  inner  product  of  two  functions  f(x)  £  L2(R )  and  g(x)  £  L2(R)  is  given  by 

/OO 

f(x)g(x)dx. 

-OO 

The  norm  of  a  function  f(x)  £  L2(R)  is  defined  as 

||/||  = 

The  convolution  of  functions  f(x)  £  L2{R)  and  g(x)  £  L2(R )  is  computed  as 

/OO 

fit)  g(x  —  t)  dt, 

-OO 

and  the  convolution  of  two  functions  f(x,  y )  £  L2(R2)  and  g(x,  y)  £  L2{R2)  equals 

/oo  r  oo 

/  fifxi  ty)  g{x  —  tx,  y  ty)dtxdty. 

-OO  j  —  00 

The  Fourier  transform  of  a  function  f{x)  £  L2{R )  is  defined  as 

/M  =  /  f(x)e~JWX  dx, 

J  —  OO 

and  the  Fourier  transform  of  a  function  f{x,y)  £  L2(R2)  is  equal  to 

/K,  uv)=  I"  /“/(*. 

J  — oo  J— oo 
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l2(Z)  and  l2(Z2)  stand  for  the  spaces  of  square-summable  discrete  signals  f(n)  and 
f{nx,ny),  respectively. 

The  ^-transform  of  a  discrete  signal  f(n)  G  l2(Z)  is  defined  as 

OO 

F(z)=  Z  f(n)z-«. 

n=— oo 

The  convolution  of  discrete  signals  f(n)  G  12{Z)  and  g(n)  G  l2(Z)  is  equal  to 

OO 

f*9{n)=  f(m)g{n-m), 

m=~  oo 

and  the  convolution  of  discrete  signals  f(nx,ny )  G  l2(Z 2)  and  g(nx,ny )  G  l2(Z2)  is  given  by 

OO  OO 

f*g(nx,ny)=  ^  £  f(mx,my)g(nx-mx,ny-my). 

TJlx  —  —  OO  ?TZy  — ■  OO 

The  Fourier  transform  of  a  discrete  signal  f(n)  G  l2(Z)  is  equal  to  the  ^-transform 
evaluated  on  the  unit  circle 

OO 

F(u)=  Z  f(n)e-iun, 

n=— oo 

and  the  Fourier  transform  of  a  discrete  signal  f(nx,  ny)  G  l2(Z2)  is  defined  as 

OO  OO 

F(u„  uy)=  Z  E  /(«..",) 

Tlx —  OO  Tly —  OO 

For  later  use,  we  define  the  following  functions: 


1.  the  unit  impulse  function 


2.  the  unit  step  function 


3.  the  rectangular  function 


for  x  =  0 
otherwise, 


for  x  >  0 
for  x  <  0, 


for  |a;|  <  | 
for  |x|  > 


4.  the  sine  function 


5.  the  unit  impulse  sequence 


where  x  G  R  and  n  G  Z. 


.  sin(7ra:) 

smc(x)  := - ,  and 


7 xx 


5(n) 


1  for  n  =  0 
0  otherwise, 
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2 


Body 

2.1  Enhancement  via  Fusion 

Existing  methods  of  mammographic  image  enhancement  can  be  divided  roughly  into  two 
categories:  (1)  methods  aimed  at  better  visualization  of  all  features  present  in  the  image 
[5,  10,  13,  14],  and  (2)  methods  that  target  specific  features  of  importance  (e.g., 
microcalcifications  [15,  16],  stellate  lesions  [9]). 

Methods  from  the  first  category  are  not  optimized  for  a  specific  type  of  cancer  and 
frequently  not  even  for  mammography.  Rather,  they  try  to  improve  the  perceptual  quality 
of  the  entire  image  and  are  often  developed  with  a  framework  more  general  than 
mammography  alone  in  mind. 

The  second  category  methods  concentrate  on  revelation  of  particular  signs  of  malignancy. 
They  can  be  very  successful  in  their  area  of  specialization;  however,  in  order  to  process 
mammogram  for  presence  of  various  features,  one  would  need  to  apply  different  algorithms 
independently  resulting  in  both  larger  number  of  images  to  be  interpreted  by  a  radiologist 
and  increased  computational  complexity  of  such  a  procedure. 

Here,  we  present  an  approach  which  overcomes  these  shortcomings  and  problematic 
limitations  via  synthesis  of  the  two  paradigms  by  means  of  image  fusion. 

The  goal  of  our  method  is  to  adapt  specific  enhancement  schemes  for  distinct 
mammographic  features,  and  then  combine  the  set  of  processed  images  into  an  enhanced 
image.  The  input  mammographic  image  is  first  processed  for  enhancement  of 
microcalcifications,  masses,  and  stellate  lesions.  From  the  resulting  enhanced  images,  the 
final  enhanced  image  is  synthesized  by  means  of  image  fusion.  Wavelet  based  image 
enhancement  and  fusion  are  merged  into  a  unified  framework,  so  that  there  is  no  need  for 
carrying  out  the  two  operations  independently  (i.e.,  computing  wavelet  decompositions, 
modifying  wavelet  coefficients  for  enhancement  of  specific  features,  reconstructing  the 
enhanced  images,  performing  wavelet  transforms  of  the  enhanced  images,  fusing  transform 
coefficients,  and  obtaining  the  final  result  by  reconstruction  from  fused  wavelet 
coefficients).  Both  enhancement  and  fusion  are  therefore  implicit  (i.e.,  performed  in  the 
wavelet  domain  only).  Figure  1  presents  a  block  scheme  of  the  overall  algorithm. 

The  algorithm  consists  of  two  major  steps:  (1)  wavelet  coefficients  are  modified  distinctly 
for  each  type  of  malignancy;  (2)  the  obtained  multiple  sets  of  wavelet  coefficients  are  fused 
into  a  single  set  from  which  the  reconstruction  is  computed.  The  devised  scheme  allows 
efficient  deployment  of  an  enhancement  strategy  appropriate  for  clinical  screening 
protocols:  enhancement  algorithm  is  first  developed  for  each  specific  type  of  feature 
independently,  and  the  results  are  then  combined  using  an  appropriate  fusion  strategy. 
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Figure  1:  Overview  of  the  algorithm. 


The  structure  of  the  algorithm  also  enables  independent  development  and  optimization  of 
enhancement  strategies  for  individual  mammographic  features  as  well  as  the  fusion  module. 


2.2  Wavelet  Transform 

In  our  previous  report,  we  modified  the  one-dimensional  discrete  dyadic  wavelet  transform 
and  used  the  obtained  representation  for  derivation  of  several  two-dimensional  transforms. 
Here,  we  focus  on  approximations  of  steerable  wavelets  and  employ  the  multiscale  spline 
derivatives  with  both  first  and  second  derivative  wavelet  decomposition.  Such  a 
representation  allows  for  both  directional  and  isotropic  processing  of  mammograms  while 
being  shift-invariant  and  aliasing-free  as  already  reported.  The  derivation  of  the  transform 
is  presented  in  Section  2.2.1.  As  is  the  case  with  the  majority  of  wavelet  transforms, 
multiscale  spline  derivatives  are  implemented  in  a  filter  bank,  so  Section  2.2.2  deals  with 
efficient  implementations  of  the  filters  in  the  filter  bank  when  the  input  image  is  mirror 
extended  to  alleviate  boundary  artifacts. 

2.2.1  Multiscale  Spline  Derivatives 

We  define  a  steerable  dyadic  wavelet  transform  of  a  function  s(x,y)  G  L2(R 2)  at  a  scale 
2m,  m  G  Z,  as  [9] 

W^ms(x,y)  =  s*^2m(x,y),  (1) 

where  i>2m(x,y)  denotes  ip2m{x,y)  rotated  by  9i,  ip2m{x,y)  =  2~2m,ip(2~mx,2~my),  ^){x,y)  is 
a  steerable  wavelet  that  can  be  steered  with  I  basis  functions,  and  9i  —  with 
i  G  {1,2,...,/}.  (For  introduction  to  steerability,  please  refer  to  our  previous  report.) 
Analogous  to  the  one-dimensional  case,  we  require  the  two-dimensional  Fourier  plane  to  be 
covered  by  the  dyadic  dilations  of  ^l(2 mu)x,  2mcoy):  there  must  exist  A3  >  0  and  B3  <  00 
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such  that 


(2) 


OO  I 

4s  <  £  £|^(2m‘‘'*,2'X)l2<  fi3 

771— —  OO  i~  1 

is  satisfied  almost  everywhere. 

If  (nonunique)  reconstructing  functions  u)  are  chosen  such  that  their  Fourier 

transforms  satisfy 

OO  I 

£  £  2mWj,)  xi(2”a.„  2mw„)  =  1,  (3) 

m— — oo  2—1 

the  function  s(x,y)  may  be  reconstructed  from  its  steerable  dyadic  wavelet  transform  by 


OO  I 

s{x,y)  =  m(x,y),  (4) 

m=-00  2—1 

where  xl m{x,y)  denotes  %2 ™{x,y)  rotated  by  0*  and  X2 ™(x,y)  =  2_2mx(2_m:r,  2 ~my). 

We  choose  wavelets  that  are  steerable  analogs  to  the  one-dimensional  derivatives  of  central 
B-spline  wavelets  [18]: 


i)(ujr,uje)  =  ( ju>rCOs(ug))d 


p+d+ 1 


(5) 


where  ujr  —  +  w2,  cog  =  a,vg(ujx,ujy),  and  d  6  {1, 2}.  These  wavelets  can  be  steered  with 

d+ 1  basis  functions. 

Wavelets  (5)  are  equal  to  d-th  order  derivatives  of  circularly  symmetric  spline  functions  in 
the  direction  of  a;- axis  (note  that  knots  for  these  splines  are  circles).  To  implement  the 
transform  efficiently,  we  approximate  the  wavelets  with  x-y  separable  wavelets 

</>(£,  y)  =  -  ^*f^A+d(y)>  (6) 

where  /3p(x)  denotes  the  central  B-spline  of  order  p. 

Based  on  the  fact  that  B-splines  tend  to  a  Gaussian  probability  density  function  as  their 
order  increases,  it  is  easy  to  see  that  both  wavelets  (5)  and  (6)  converge  to  the  same 
functions  (i.e.,  d-th  order  derivatives  of  the  normalized  Gaussian  in  the  direction  of  x-axis) 
as  p  — >■  oo.  In  order  to  steer  wavelets  tp(x,  y)  given  by  (6)  (note  that  steering  will  be  only 
approximate,  since  these  wavelets  are  not  steerable),  we  need  to  find  basis  functions  that 
will  approximately  steer  ip(x,y).  To  accomplish  this,  we  take  advantage  of  the  property  of 
circularly  symmetric  functions  that  rotations  of  their  directional  derivatives  are  equal  to 
directional  derivatives  in  rotated  directions: 

j  ddgc(x,y)}  ddgc(x,y) 
e°\  dnd  }~  dfii  ’ 
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where  TZg0  stands  for  rotation  by  $o,  =  nV Qc(x,y),  Qc{x,y)  is  a  circularly  symmetric 

function,  and  ng0  denotes  vector  n  =  (cos  9,  sin  9)  rotated  by  90. 

Let  us  choose 

e(x,y)  =  /3p+d(x)/3p+d(y), 

which  is  approximately  circularly  symmetric  function  for  higher  order  splines.  A  rotation  of 
ij){x,y)  =  from  Equation  (6)  by  90  can  therefore  be  approximated  by 


il>e°(x,y)  ~ 


y) 

dnd 


d-ii  %+djx )  <PPp+d{y) 

x  y  dxd_i  dy{ 


(7) 


where  n  =  (cos#0,  sin^0)  =  (nx,ny).  (Note  that  in  case  of  Gaussian,  which  is  both  x-y 
separable  and  circularly  symmetric,  Equation  (7)  becomes  exact.) 

To  derive  an  algorithm  for  the  fast  computation  of  the  transform,  we  introduce  two 
smoothing  functions  such  that 


oo  I 

fax,  cjy)  fax,  s)=EE  n2mu>x,  2mujy)  xi(2mo;x,  2muy).  (8) 

m~ 0  i= 1 

Using  a  set  of  basis  functions  (7)  that  approximately  steer  wavelets  (6),  we  want  to 
construct  a  transform  such  that  Equations  (1)  through  (4)and  (8)  will  be  valid  (superscript 
i  must  be  viewed  now  as  an  index,  rather  than  rotation  by  9t). 

Let  F(uj)  be  a  digital  filter  frequency  response  and  let  us  denote 


Fs(cu)  =  ei“sF(uj) 


with  s  being  a  filter  dependent  sampling  shift  needed  to  obtain  finite  impulse  response 
(FIR)  filters. 

In  frequency  domain,  we  can  express  basis  functions  from  (7)  as 

ipt+l(u)x,ujy)  =  Gi~/(cox)Gl_s(Ljy)Pp+i(ux)pp+dfay),  i  €  {0, 1, 2},  (9) 

where  Gd{u> )  is  given  by 

(?>(«)  =  e*"(2j  sin  (1))'.  (10) 

where  d  is  the  order  of  the  derivative,  d  £  {1,2},  the  sampling  shift  for  filter  (10)  is 
s  =  and  G°(cj)  =  1. 

Since  we  are  interested  in  the  first  and  second  derivative  wavelets,  we  impose 

<p(cjx,  U>y)  =  (fi(u>x,  U>y)  =  /3p(u>x)ftp(u)y) 
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and  choose 


X  =  Ks(^x)T(u)y)(3p(u)x)Pp—2{u>y), 

X  (WIl  W!l)  =  ■^s(Li>x)Ks(u)y)Pp-l(U>X)^p-l(u)y), 

X  (wx,0}y)  =  T(UX)Ks(c0y)Pp-2(U>x)flp(t0y), 


(11) 

(12) 

(13) 


where 


with  the  sampling  shift  for  Kd(u>)  being  the  same  as  the  one  for  Gd(uj),  and 

T(w)  =  \H(u)\2 


(14) 


(15) 


H(w)  =  «*■»  (cos  (|))”+1  ■ 

the  sampling  shift  for  H(u>)  being  s=  (?H~1)^nod2, 

Using  the  relation  derived  in  the  previous  report 

Pp(  2w)  =  H-a(u)Pp(w) 


(16) 


(17) 


together  with  Equations  (9)  and  (11)  through  (13)  with  Equation  (8)  results  in 

G2{ux)K2(u>x)T(u>y)  +  G1  {ux)Kl  {ux)Gl  {u)v)Kl  (uy)  +  T(uJx)G2(uy)K2(uJy)+ 

+\H(cox)H(uy)\2  =  1. 


Next,  we  derive  a  filter  bank  implementation  of  the  transform.  Assuming  a  bandlimited 
input  signal  s(u>x,u>y)  =  0  for  |o>x|  >  7r  or  \u)y\  >  it  and  using  Shannon’s  sampling  theorem  in 
two  dimensions  [19]  with  Equation  (1)  and  basis  functions  from  Equation  (9),  we  can  write 


roc  roc  °o^ 

is(x,  y)=  /  Y  X  s(^>  *v)  sinc(^  -  »*)  sinc(ty  -  *„)• 

J—OOJ—OO;  _ _ •  _ _ 


W\ 


00  00  ix=~oo  iy——oc 


■  X  9-st(,mx)PP+i(x  -tx-  mx)  Yj  9-s{my)Pp+d-i{y  ~ty-  my)  dtx  dty , 

mx=— oo  my——oc 


where  i  6  {0, 1, 2}  as  in  Equation  (9). 

We  approximate  sine  functions  with  r- order  cardinal  splines,  then  use  the  relation  between 
cardinal  and  B-splines 

OO 

Vr(x)=  X  K'^Prix-i) 

i=- oo 
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and  get 

T{W2ms(x,y)\  }  ~  S{ux,u)y)  B~x(u)x)  B~l{uiy)  Bp+r+i+1(u)x)- 

I  x=nXyy=Tiy 

m— 1 

••Bp+rW- w(u>y)  Gt;( -Tu.)  G1,(2"V,)  n  Hp-ti(rUx)Hitd-t(2’'wv),  (18) 

71=0 

where  B~1{oj)  denotes  the  Fourier  transform  of  the  direct  B-spline  filter  of  order  p.  Table  1 
shows  the  ^-transforms  of  direct  B-spline  filters  for  the  first  ten  orders. 

Using  Equation  (18)  with  an  approximation  Bp+r+i+i(u> )  ~  Bp+r(u>)Bl (a;),  we  can  obtain  a 
filter  bank  implementation  of  the  transform  decomposition.  The  reconstruction  part 
follows  from  Equations  (8),  (9),  and  (11)  through  (13).  Figure  2  shows  a  filter  bank 
implementation  of  the  transform.  Noninteger  shifts  at  scale  1  are  rounded  to  the  nearest 
integer.  Tables  2  through  5  list  impulse  responses  of  the  filters  used  in  the  filter  bank  for 
P  €  {0, 1, 2}. 

The  derived  transform  enables  both  second  derivative  directional  analysis  and  Laplacian  of 
Gaussian  approximations  across  dyadic  scales  (the  latter  can  be  obtained  through 
summation  of  the  outputs  from  blocks  (?2(2ma>)  applied  along  x  and  y  axis).  Furthermore, 
addition  of  a  block  Gx_s(2mu)  at  each  level  of  decomposition  allows  first  derivative 


16 


B>x)  B'rCCOy)  - 

-  Bp+r+1(cox)Bp+r+1(coy) 

(a) 

Bp’+r+l(®x)  BpVr+1(COy)  - 

-  Br((ox)  Br(coy) 

(b) 


(c)  (d) 


Figure  2:  Filter  bank  implementation  of  multiscale  spline  derivatives  for  m  G  [0,  M  —  1]:  (a) 
Prefiltering,  (b)  postfiltering,  (c)  decomposition  module,  and  (d)  reconstruction  module. 


Table  2:  Impulse  responses  px(n)  and  g2(n). 

f(n) 


n 


-1 

0 

1 


(/(n) 


1 

-1 


Table  3:  Impulse  res 


e  res 

ponses 

h(n),  k 

1(n),  k 2 

(n),  an 

n 

h(n) 

fc1(n) 

k2(n) 

t(n) 

-i 

ia 

0.25 

0 

EH 

-0.25 

-0.25 

0.5 

1 

0.25 

0.25 

Table  4:  Impulse  responses  h(n ),  kl(n),  k2(n),  and  t(n )  for  p  =  1. 


n 

h(n) 

A;1(n) 

k2(n) 

t(n) 

-2 

-1 

0.25 

-0.0625 

-0.0625 

0.0625 

0.25 

0 

0.5 

-0.3125 

-0.375 

0.375 

1 

0.25 

0.3125 

-0.0625 

0.25 

2 

0.0625 

0.0625 
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Table  5:  Impulse  responses  h(n),  Ad(n),  k2(n),  and  t(n)  for  p=2. 


n 

h(n ) 

kl(n ) 

k2(n) 

t{n) 

-3 

-2 

0.125 

-0.015625 

-0.015625 

0.015625 

0.09375 

-1 

0.375 

-0.109375 

-0.125 

0.234375 

0 

0.375 

-0.34375 

-0.46875 

0.3125 

1 

0.125 

0.34375 

-0.125 

0.234375 

2 

0.109375 

-0.015625 

0.09375 

3 

0.015625 

0.015625 

(a)  (b) 


Figure  3:  Spline  derivatives  in  the  rr-axis  direction,  (a)  Wavelet  equal  to  the  first  derivative 
of  a  quartic  spline,  (b)  Wavelet  equal  to  the  second  derivative  of  a  quintic  spline. 


directional  analysis  as  well.  Figure  3  shows  first  and  second  derivative  wavelets  obtained  as 
linear  combinations  of  cubic  B-splines. 

2.2.2  Filter  Implementations 

Since  all  two-dimensional  filters  used  in  the  filter  bank  implementation  of  the  transforms 
are  x-y  separable,  only  one-dimensional  filters  need  to  be  implemented.  We  describe  the 
implementation  of  finite  impulse  response  (FIR)  filters  first  and  then  treat  prefiltering  and 
postfiltering  infinite  impulse  response  (HR)  filters  from  Figure  2(a)  and  (b). 

Let  us  refer  to  filters  applied  at  scale  2m  as  filters  at  level  rn+1,  and  let  filters  at  level  1 
(Equations  (10),  (14),  (15),  and  (16)  be  called  “original  filters,”  to  distinguish  them  from 
their  upsampled  versions.  Let  us  split  the  input  to  the  filter  bank  from  Figure  2  into  image 
matrix  rows  and  columns,  each  corresponding  to  a  real  signal  s(n)  6  12{Z),  n  €  [0,  N  —  1]. 
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Depending  on  the  length  of  each  filter  impulse  response,  filtering  an  input  signal  may  be 
computed  either  by  multiplying  the  discrete  Fourier  transforms  of  the  two  sequences  or  by 
circularly  convolving  s(n)  with  a  filter’s  impulse  response.  Using  circular  rather  than  linear 
convolution,  as  is  customary  in  image  processing,  can  lead  to  boundary  artifacts  caused  by 
abrupt  changes  in  the  periodically  extended  signal.  A  common  remedy  for  such  a  problem 
is  realized  by  constructing  a  mirror  extended  signal  [12] 


Sme(^)  — 


s(—n  —  1)  if  n  G  [-N,  —1] 
s(n)  if  n  G  [0,  N  —  1], 


(19) 


where  we  chose  the  signal  sme(n )  to  be  supported  in  [— N,  N  —  1]. 

Let  us  classify  symmetric/antisymmetric  real  even-length  signals  into  four  types  [20]: 


Type  I  f(n)  =  f(-n), 

Type  II  f(n )  =  f(-n  -  1), 

Type  III  /(n)  =  -f(-n), 

Type  IV  f(n)  =  -f(-n  -  1), 

where  n  G  [— N,  N  —  1].  Note  that  for  Type  I  signals  the  values  at  /( 0)  and  f(—N )  are 
unique,  and  that  for  Type  III  signals  the  values  at  /( 0)  and  f(—N )  are  equal  to  zero. 

(This  is  important  for  storage  requirements:  for  signals  of  Type  II  or  Type  IV,  N  samples 
need  to  be  saved,  while  Type  I  and  Type  III  signals  require  N+l  and  N— 1  sample 
representations,  respectively.) 

Using  properties  of  the  Fourier  transform,  it  is  easy  to  show  that  the  convolution  of 
symmetric/antisymmetric  real  signals  results  in  a  symmetric/antisymmetric  real  signal.  If 
a  symmetric/antisymmetric  real  signal  has  an  even  length,  then  there  always  exists  an 
integer  shift  such  that  the  shifted  signal  belongs  to  one  of  the  above  types. 

Now,  we  are  ready  to  examine  the  filter  bank  implementation  of  the  wavelet  transform 
from  Figure  2  with  filters  given  by  Equations  (10),  (14),  (15),  and  (16)  driven  by  mirrored 
signals  of  the  form  sme(n )  from  Equation  (19)  at  the  input.  Let  the  number  of  levels  M  be 
restricted  by 

M  <  1+log,  .iV~_11,  (20) 

L'max  -L 

where  Lmax  is  the  length  of  the  longest  original  FIR  filter  impulse  response. 

Each  FIR  filter  block  in  the  filter  bank  consists  of  a  filter  and  a  circular  shift  operator. 
Equation  (20)  guarantees  that  the  length  of  the  filter  impulse  response  does  not  exceed  the 
length  of  the  signal  at  any  block. 
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Since  our  mirror  extended  input  row  or  column  sme(n)  is  of  Type  II  and  noninteger  shifts 
at  level  1  are  rounded  to  the  nearest  integer,  it  follows  that  a  processed  one-dimensional 
signal  at  any  point  in  the  filter  bank  belongs  to  one  of  the  types  defined  above.  This  means 
that  filtering  a  signal  of  length  2 N  can  be  reduced  to  filtering  a  signal  of  approximately  one 
half  of  its  length. 

Implementation  is  particularly  simple  for  FIR  filters  designed  with  2  and  p  odd.  Filters 
are  of  Type  I  in  this  case,  so  their  output  will  be  of  Type  II.  An  FIR  filter  block  from  the 
filter  bank  shown  in  Figure  2  can  therefore  be  implemented  by 

L  — 1 

Fs,mu(n)  =  /(0)it//(n)  +  ^  /00[u//(n  ~  2m*)  +  un(n  +  2m*)]>  n  G  [0,  N  —  1],  (21) 

i—  1 


where 

f  u(-n-l)  ifne[-y,-l] 

Un{n)  =  <  u\n)  if  n  G  [0,  N  —  1]  (22) 

[  u(2iV  —  n  —  1)  if  n  G  [N,  &-], 

u(n)  is  an  input  signal  to  a  block,  f(n)  is  an  impulse  response  of  G2( 2ma>),  K2(2mu), 
T(2mu),  or  H{2muj)  with  p  odd,  L  is  the  length  of  the  filter,  and  N  is  the  length  of  an 
input  signal  s(n)  to  the  filter  bank.  Implementation  of  filters  bp(n)  used  for  prefiltering  and 
postfiltering  (Figure  2(a)  and  (b))  represents  a  special  case  of  Equation  (21)  with  m  =  0. 

A  filter  bank  with  the  above  implementation  of  blocks  and  signal  s(n )  at  the  input  yields 
equivalent  results  as  circular  convolution  of  input  sme(n )  as  defined  by  Equation  (19).  In 
addition  to  requiring  one  half  the  amount  of  memory,  the  computational  savings  over  a 
circular  convolution  implementation  of  blocks  are,  depending  on  the  original  filter  length, 
three  to  four  times  fewer  multiplications  and  one  half  as  many  additions. 

A  similar  approach  is  used  for  other  filters.  The  problem  becomes  slightly  more  involved  in 
this  case,  because  the  filters  change  type  from  first  to  subsequent  levels,  and  the  signal 
component  type  can  be  altered  by  a  filter  block  as  well.  As  a  consequence,  an 
implementation  of  blocks  that  use  distinct  original  filters  may  not  be  the  same,  and  the 
implementation  of  blocks  at  level  1  may  differ  from  the  implementation  of  blocks  at  other 
levels  of  analysis. 

The  decomposition  blocks  at  level  1  can  be  implemented  by 


and 


i-1 

G-s,ou(n)  =  Y,  9(i)[un{n  -  *  -  1)  -  un(n  +  *)],  n  e  [1,  N  -  1] 
i= 0 


L- 1 
2  1 

H-Sfiu(n)  =  ^/i(i)[u//(n-i-l)  +  tt//(n  +  i)],  n€[0,IV], 

z=0 
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for  p  even,  where  ujj(/ )  is  defined  by  (22),  g(n)  and  h(n)  are  impulse  responses  of  the 
filters  computed  from  (10)  and  (16),  respectively,  and  L  is  the  length  of  the  corresponding 
impulse  response. 

The  output  from  a  block  GLs(w)  at  level  1  is  of  Type  III,  while  the  output  from  H-S(u>)  at 
the  same  level  is  of  Type  I. 

The  decomposition  blocks  at  subsequent  levels  m  G  [1,M— 1]  can  be  implemented  by 


G-s,mu(n )  =  Z  g(i)[Mn  -  2m(*  +  s))  -  Mn  +  2m(*  +  «))].  n  G  [1,  N  -  1], 

i= 0 


for  p  even, 


u-i 

G-S,mu(,n)  =  Z  9{i)[un(n  ~  2 m(i  +  s ))  -  un(n  +  2 m(i  +  «))],  n  G  [0,  N  -  1], 
1=0 


for  p  odd, 


L  — 1 

F_s>TOu(n)  =  /(0)it7(n)  +  Z  /(*)[«/(”  ~  2m«)  +  u/(n  +  2mi)],  nG[0,iV],  (23) 

2=1 


with  /(n)  =  g(n)  for  d= 2  and  p  even, 


H-S,mu(n)  =  Z  M0[u/(n  “  2m(*  +  s))  +  «/(n  +  2m(i  +  s))],  n  G  [0,  iV],  (24) 
2=0 


for  p  even,  where 

f  u(-n )  if  n  G  [— y,  -1] 

U/(n)  =  <  u(n)  if  n  G  [0,  N]  (25) 

k  u(2iV  —  n)  if  tig  [N+l,^]. 

The  outputs  from  blocks  G-s( 2mu)  are  of  Type  III  for  d=l  and  p  even,  of  Type  IV  for 
d—  1  and  p  odd,  and  of  Type  I  for  d=2  and  p  even,  whereas  the  outputs  from  H^s(2mu>) 
are  of  Type  I  for  p  even. 

Next,  the  reconstruction  blocks  at  level  1  can  be  implemented  by 


L 

2 

K]fiu{n)  =  Z  k{i)[uni{n  -  i  +  1)  -  uin(n  +  t)],  n  G  [0,  N  -  1] 
2=1 


and 

L 

2 

HSt0u(n)  =  Z  M*)[M/(n  —  i  +  1)  +  ui(n  + 1)],  n  G  [0,  N  —  1], 

i=l 
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for  p  even,  where 


—u(—n) 

0 


uniin)  =  { 


u(n) 

0 


[  —u(2N  —  n) 


if  ne  [— f,-l] 

if  n  =  0 

if  n  G  [1,  iV  —  1] 
if  n  =  N 

ifn€  [N  +  l,*!f\, 


(26) 


ui(n)  is  as  defined  by  (25)  and  k(n)  is  an  impulse  response  of  the  filter  from  (14).  Note 
that  both  outputs  from  blocks  K\  (a;)  and  Hs(oo)  are  of  Type  II. 

The  reconstruction  blocks  at  subsequent  levels  can  be  implemented  by 


2  x 

Ks,mu(n )  =  H  *(*  +  !)[«///(«  “  2m(?  +  a))  -  uIH(n  +  2 m{i  +  a))],  n  G  [0 ,N], 

2—0 

for  d=  1  and  p  even,  (23)  with  /(n)  =  A:(n)  for  d— 2  and  p  even, 


2  1 


Kdsmu{n)  =  ^2  k(i  +  l)[uIV(n  -  2 m(i  +  a))  -  «jv(n  +  2m(i  +  a))], 
2—0 


n  G  [0,  iV  —  1], 


for  d—1  and  p  odd, 


for  p  even,  where  um(l)  is  given  by  (26), 


uiv(n ) 


—u(—n  —  1)  if  n  G  [—  y,—  1] 

<  u(n )  if  n  G  [0,  N  —  1] 

—u(2N  —  n  —  1)  if  n  G  [N,  ^], 


and  H-StTnu(n)  is  specified  by  Equation  (24).  We  observe  that  the  outputs  from  blocks 
Kf(2muj )  and  Hs(2mco),  m  G  [1,  M  —  1],  for  p  even  are  of  Type  I. 

When  we  compare  the  above  implementation  of  blocks  to  circular  convolution  driven  by  a 
mirrored  signal  sme(n )  at  the  input,  we  observe  that  approximately  twofold  less  memory 
space,  three  to  four  times  fewer  multiplications  and  one  half  as  many  additions  are 
required.  (For  Type  I  signals  an  additional  sample  has  to  be  stored  because  two  values  are 
without  a  pair) . 

The  implementation  presented  in  this  section  performs  all  operations  in  the  spatial  domain; 
however,  one  could  also  implement  the  structures  shown  in  Figure  2  with  an  input  signal 
sme(n)  in  the  frequency  domain.  For  short  filter  impulse  responses,  such  as  those  given  in 
Tables  3,  4  and  5,  the  spatial  implementation  described  in  this  section  is  certainly  more 
efficient.  For  long  filter  impulse  responses,  however,  filtering  is  faster  if  implemented  in  the 
frequency  domain.  Additional  details  on  alternative  FIR  filter  implementation  strategies 
can  be  found  in  [21]. 
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Implementation  of  HR  filters  b~ 1(n)  used  for  prefiltering  and  postfiltering  is  a  bit  more 
involved  than  the  one  of  their  FIR  counterparts.  Fortunately,  the  number  of  different  cases 
is  much  smaller  here:  possible  input  to  b~ *(n)  in  the  filter  bank  from  Figure  2  is  either  of 
Type  II  or  of  Type  I  (symmetry  types  for  HR  filters  slightly  differ  from  those  defined  for 
FIR  filters:  here,  mirror  extended  signals  are  periodically  repeated,  so  that  they  stretch 
from  —oo  to  oo).  We  use  ideas  and  a  few  results  from  [22]. 

Let  us  first  take  a  closer  look  at  the  system  function  B~1(z)  with  p£  {2,3}.  This  function 
can  be  written  as  a  cascade  of  terms 


.  .  i  —  a 

Z  —  Z~X  (1  — az_1)(l  —  az)’ 

which  can  be  expressed  in  a  parallel  form  as 

£(*)  =  (  1  t  - - l)  , 

1  —  or  V 1  —  az  1  1  —  az  / 

where  a  and  L  are  poles  of  the  causal  and  the  anticausal  filter,  respectively. 
The  impulse  response  of  this  term  can  be  written  as 

e(n)  =  T~^alnl- 

1  —  or 


(27) 


(28) 


We  choose  to  implement  E(z)  in  a  cascade  form  and  therefore  extract  the  difference 
equations  from  Equation  (27): 


c+(n)  =  u(n)  +  a  c+(n  —  1)  n  =  1, 2, ... ,  N—l , 


(29) 


and 

c(n)  =  a  (c(n  +  1)  —  c+(n))  n  =  N— 2,  N— 3, . . . ,  0,  (30) 

where  u(n)  denotes  the  input  to  the  block,  c+(n)  is  the  output  from  the  causal  part,  and 
c(n)  stands  for  the  output  from  the  block. 

To  solve  Equations  (29)  and  (30)  we  need  boundary  conditions  c+(0)  and  c(N—  1).  We 
derive 

0  N- 1  i+ 1  ,  fJiN-i  io 

c+(°)  =  E  =  w(°)  +  E  1  _  Ton  -  u(0)  +  £a<+1«(*)» 

i— 2—0  1  ^  -•—I n 


(31) 


i=0 


and,  using  parallel  form  (28) 


<N-i)  =  r-^r2(c*(N-i)+Y: 


»--l  aK~i  +  Q»+l+i 


1  —  a2 

—a 
1  -  a5 


l-a2N 
(c+(N- 1)+  E  aN~lu(i))i 


~u 


(0) 


i= 0 
N- 1 


(32) 


2— TV -1-20 
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where 


unp(n)  = 


uu(n  mod  ( 2N ))  if  n  >  0 

un(—(n  +  1)  mod  ( 2N ))  if  n  <  0, 


un(n)  = 


u(n)  if  n  E  [0,  N  —  1] 

u(2N  —  n  —  1)  if  n  €  [N,  2N  —  1], 


N  is  the  length  of  an  input  signal  to  the  filter  bank,  and  io  <  N— 1  is  selected  such  that  a10 
falls  below  a  predefined  precision  threshold. 

For  orders  p  greater  than  three,  we  implement  B~1(z)  as  a  cascade  of  terms  E(z)  with 
different  a’s.  Note  that  the  output  from  block  E(z)  is  always  of  the  same  type  as  the  input 
to  it. 


2.3  Enhancement  of  Mammographic  Features 

2.3.1  Microcalcifications 

Microcalcifications  appear  on  mammograms  in  approximately  half  of  breast  cancer  cases. 
The  assessment  of  shape,  number,  and  distribution  of  microcalcifications  is  important  for  a 
radiologist  to  reach  the  correct  diagnosis.  Microcalcifications  are  smaller  than  1  mm  in  size 
and  can  be  difficult  to  locate  when  they  are  superimposed  on  dense  breast  tissue. 

Several  techniques  have  been  developed  to  improve  the  visibility  of  microcalcifications 
[15,  16,  23,  24].  The  approach  devised  by  Strickland  and  Hahn  [15]  is  particularly  well 
suited  for  our  framework:  they  used  an  undecimated  wavelet  transform  to  approximate 
second  derivatives  of  a  Gaussian  probability  density  function  for  a  multiscale  matched 
filtering  for  presence  of  microcalcifications. 

Strickland  and  Hahn  based  their  method  on  the  observation  that  the  average 
microcalcification  can  be  modeled  by  a  circularly  symmetric  Gaussian  function.  Using  a 
combination  of  a  separable  Markov  process  with  autocorrelation  rnn  =  o2ne  aVw+\f\  and  a 
nonseparable  Markov  process  with  autocorrelation  rnn  =  cr2e~a^k2+l2  to  represent 
mammogram  texture,  they  obtained  the  separable  matched  filter 

MSep(UX,UJy)  =  M(u>X)M(u}y),  (33) 


where 

M(lo) 

and  the  nonseparable  matched  filter 


OJKU1 

- - - g  2 

V2aan 


(34) 
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In  order  to  deal  with  different  sizes  of  microcalcifications,  one  must  vary  a  of  matched 
filters  (33)  and  (34)  appropriately.  In  Strickland  and  Hahn’s  scheme,  a  wavelet 
decomposition  was  chosen  to  approximate  the  matched  filters  across  the  desired  scale 
range.  Considering  the  100/xm  resolution  of  the  Nijmegen  database,  the  wavelet  transform 
was  computed  over  the  first  4  octaves.  For  a  denser  sampling  of  scale,  voices  were  inserted 
at  octaves  “2.5”  and  “3.5.”  The  wavelet  analysis  stage  acted  as  a  bank  of  matched  filters; 
wavelet  coefficients  at  locations  indicating  microcalcifications  were  multiplied  by  a  gain 
factor,  and  then  the  inverse  wavelet  transform  was  applied  to  the  modified  coefficients. 

In  our  approach,  microcalcifications  are  modeled  by  central  B-splines.  Using  the  relation 
between  the  standard  deviation  of  a  Gaussian  function  and  the  order  of  B-splines  that 
approximate  it  a  =  [25],  the  assumption  that  a  Gaussian  object  is  visible 

approximately  over  ±a  pixels  [15],  and  the  fact  that  the  mammograms  in  the  University  of 
Florida  database  were  digitized  at  116/zm  resolution,  four  levels  of  the  transform  from 
Section  2.2  with,  for  example,  p= 3  are  needed  to  encompass  different  sizes  of 
microcalcifications.  The  wavelet  decomposition  including  voices  at  scales  3  and  6 
(corresponding  to  Strickland  and  Hahn’s  octaves  “2.5”  and  “3.5”)  can  be  obtained  by 
deriving  a  counterpart  to  Equation  (18)  for  the  two  scales. 

0p(3w)  can  be  related  to  (3p(lu )  by  expressing  0P(3 to)  as  (cf.  Proposition  1  of  [26]) 


0p(  3w)  = 


1 

3P+1 


/'sm(f)\',+1  /  sin  (f )  \ 

U"(f)j  l  t  J 


Using  we  get 


Pp{  3w)  =  V(w)0p(w), 


where 

V(u)  =  (i^  +  l  +  e*'))’’  . 

Filter  V(uj)  can  be  implemented  as  a  moving  sum  with  2(p+l)  additions  per  sample  and  a 
multiplicative  factor  applied  to  the  wavelet  coefficients  [26] . 

Next,  0P( 6o»)  is  expressed  by  means  of  Equation  (17): 

4,(6o;)  =  H-s(3co)/3p(3uj) 


with  odd  orders  p  used  for  the  sampling  shift  s  to  be  zero. 
Now,  we  can  write 


?{Wls{x,y) I  }  ~  S(ux,u)y)Br  1  (ux)  Br  1(uy)  Bp+r+i+1(ujx)- 

1  x=nx,y=ny 

'Bp+T+d-i+i(uy)  Gf '( 3ux)  Gi_,(3uy)VT+iMV+d-iM 


(35) 
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and 


T{W£s(x,y)  }  ~  S{ujx,uy)  Br  ‘(u,)  Br  1(uy)  Bp+r+i+1(ujx)- 

x=nXiy=ny 

'Bp+r+d-i+i(vy)  GtA  6w.)  GU(6w„)J/^‘(3wa!)JK:''-i(H)',P+‘(^)^’+,‘"i("!.)  (36) 


with  notation  being  the  same  as  in  Equation  (18),  and  superscript  p  in  Hp_s{ u>)  denoting 
the  order  p  in  Equation  (16). 

Wavelet  coefficients  obtained  via  Equations  (35)  and  (36)  are  not  used  for 
reconstruction — the  inverse  transform  is  carried  out  as  given  in  Section  2.2. 

The  decomposition  described  by  Equations  (18),  (35),  and  (36)  with  additional  filtering  by 
G2(lu) )  at  each  scale  l  €  {1,2,3, 4, 6, 8}  enables  approximations  to  the  second  derivatives  of 
Gaussian  along  both  x  and  y  directions  and  to  Laplacian  of  Gaussian  across  distinct  scales 
employed  by  Strickland  and  Hahn  [15]  (cf.  Equations  (33)  and  (34)).  We  proceed  in  a 
similar  fashion  as  therein:  the  two  outputs  per  scale  are  thresholded  independently,  all 
binary  results  are  then  combined,  a  circular  region  centered  at  detected  pixel  locations  are 
next  multiplied  by  a  gain,  and,  finally,  the  reconstruction  process  uses  modified  transform 
coefficients. 


2.3.2  Circumscribed  Masses 


Almost  half  of  missed  cancers  appear  on  mammograms  as  masses.  Perception  is  a  problem 
particularly  for  patients  with  dense  fibroglandular  patterns.  The  detection  of  masses  can 
be  especially  difficult  because  of  their  small  size  and  subtle  contrast  compared  with  normal 
breast  structures. 

Fan  and  Laine  [10]  developed  a  discrete  dyadic  wavelet  transform  based  algorithm  suitable 
for  enhancement  of  masses.  They  constructed  an  approximation  to  Laplacian  of  Gaussian 
across  dyadic  scales  for  an  isotropic  input  to  a  piecewise  linear  enhancement  function. 
Approximation  to  Laplacian  of  Gaussian  across  dyadic  scales  is  easy  to  obtain  using 
multiscale  spline  derivatives  derived  in  Section  2.2:  Equation  (18)  with  i  =  0  and  i  =  2 
approximates  the  second  derivative  of  a  Gaussian  function  along  directions  of  x  and  y  axis, 
respectively  (the  corresponding  branches  in  Figure  2(c)  are  the  first  and  third  from  the 
top).  The  appropriate  transform  coefficient  at  each  dyadic  scales  are  therefore  added  and 
their  sum  input  to  the  piecewise  linear  function 


C(x)  =  { 


x  -  (K  -  1  )T 
Kx 

x  +  (K-  1  )T 


if  x  <  —T 
if  |z|  <  T 
if  x  >  T. 


(37) 


used  at  each  level  m+1  of  the  transform  separately.  Due  to  the  expected  size  of  masses, 
levels  greater  than  4  are  enhanced  more  aggressively. 


26 


25 


Figure  4:  Enhancement  function  (Equation  (37)  with  K  =  20  and  T  =  1). 

Figure  4  shows  the  enhancement  function  from  Equation  (37)  for  parameter  values  K  =  20 
and  T  =  1. 

The  multiplicative  factor  obtained  as  the  ratio  between  the  output  and  input  of  the 
enhancement  function  is  next  applied  to  the  original  wavelet  coefficients  [10],  and  then  the 
reconstruction  (Figure  2(b)  and  (d))  is  carried  out. 

Figure  5  shows  the  cranio-caudal  view  of  a  patient  with  bloody  nipple  discharge.  On  the 
enhanced  image  cropped  to  the  area  of  interest,  irregular  anterior  borders  of  a  mass  are 
better  seen. 

2.3.3  Stellate  Lesions 

It  is  important  for  radiologists  to  identify  stellate  lesions  since  their  presence  is  a  serious 
indicator  of  malignancy.  Stellate  lesions  vary  in  size  and  subtlety  and,  in  addition,  do  not 
have  a  clear  boundary,  making  them  difficult  to  detect. 

In  the  development  of  our  algorithm,  we  follow  an  observation  made  by  Kegelmeyer  et  al. 
about  the  distortion  of  edge  orientation  distribution  induced  by  a  stellate  lesion  [27]. 
Normal  mammograms  show  a  roughly  radial  pattern  with  structure  radiating  from  the 
nipple  to  the  chest  wall.  A  stellate  lesion  not  only  changes  this  pattern,  but  also  creates 
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Figure  5:  Contrast  enhancement  of  the  cranio-caudal  of  a  patient  with  bloody  nipple  dis¬ 
charge.  (a)  The  original  mammogram  with  area  of  interest  delineated,  (b)  Unprocessed 
extracted  area,  (c)  The  enhanced  area  improves  the  visualization  of  the  mass. 
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another  center  from  which  rays  radiate.  Directional  analysis  using  the  Sobel  edge  operator 
was  employed  for  assessment  of  local  orientations  [27]. 

Wavelet  transform  from  Section  2.2  enables  directional  analysis  as  well.  By  adding 
additional  filter  G1_s(2mu))  to  each  scale  of  decomposition  from  Figure  2,  approximations  to 
both  first  and  second  steerable  derivatives  of  a  Gaussian  are  available.  A  multiscale 
derivative-pair  quadratic  feature  detector  is  computed  by  finding  the  maximum  of  the  local 
oriented  energy  with  respect  to  angle  9 

Ee2m  (x,y)  =  V(VFl^s(rr)2/))2  +  (W2^S(x,y))2,  (38) 

where  W\e2ms{x,  y)  and  W2e2ms(x,  y)  denote  wavelet  decompositions  using  first  (Equation 
(6)  with  d=  1)  and  second  (Equation  (6)  with  d=  2)  derivative  wavelet,  respectively, 
steered  to  angle  6.  The  angle  that  maximizes  the  local  oriented  energy  (38)  represents 
orientation  at  pixel  location  (x,y). 

Similarly  to  the  method  from  Section  2.3.1,  processing  is  carried  out  within  windows  with 
scale  dependent  sizes:  1-norm  of  differences  between  the  local  and  average  orientations  is 
computed  in  the  window  and  used  as  a  measure  of  orientation  nonuniformity.  Soft 
thresholding  as  a  function  of  the  orientation  nonuniformity  measure  is  next  applied  to  the 
transform  coefficients  at  each  dyadic  scale  independently  [9].  The  altered  coefficients  are 
then  included  for  reconstruction. 

Figure  6  shows  the  oblique  view  with  a  mass  visible  in  the  mid-posterior  breast.  The 
enhanced  image  demonstrates  the  irregularity  and  spiculation  of  the  mass. 

2.4  Fusion  of  Enhanced  Features 

Enhanced  features  from  Section  2.3  are  fused  into  the  final  enhanced  image.  Image  fusion 
using  redundant  steerable  wavelet  representations  was  treated  in  detail  in  our  previous 
report,  so  let  us  only  point  out  that  the  fusion  used  here  is  implicit  rather  than  explicit: 
wavelet  transform  coefficients  are  first  modified  for  enhancement  of  microcalcifications, 
circumscribed  masses,  and  stellate  lesions,  and  then  the  new  coefficients  are  obtained  by 
fusion  before  the  reconstruction  is  accomplished. 

Note  also  that  it  is  possible  to  put  different  weights  on  features,  and  exclude  certain 
features  from  the  final  result. 
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Figure  6:  Contrast  enhancement  of  the  oblique  view,  (a)  The  original  mammogram  with 
area  of  interest  delineated,  (b)  Unprocessed  extracted  area,  (c)  The  enhanced  area  improves 
the  visibility  of  the  mass. 


3  Conclusions 


During  our  second  year,  we  developed  a  method  that  improves  the  visibility  of  specific 
mammographic  features  by  local  enhancement  and  fusion  of  the  enhanced  areas.  The 
described  method  incorporates  a  variety  of  properties  of  mammographic  image 
enhancement  methods  tailored  to  specific  signs  of  malignancy  into  a  unified  computational 
framework. 

We  derived  a  wavelet  transform  with  analysis  stage  enabling  approximations  to  directional 
first  and  second  derivatives  of  a  Gaussian  function  and  to  Laplacian  of  Gaussian  across 
distinct  scales.  Such  a  decomposition  is  suitable  for  both  anisotropic  and  isotropic 
multiscale  analysis,  and,  in  addition,  provides  an  adaptable  framework  for  incorporation  of 
a  variety  of  mammographic  processing  methods. 

The  derived  transform  has  also  proved  flexible  enough  for  enhancement  and  fusion  of 
individual  types  of  mammographic  features.  Separate  enhancement  algorithms  have  been 
developed  for  microcalcifications,  circumscribed  masses,  and  stellate  lesions,  and  fusion  of 
the  modified  transform  coefficients  performed  before  the  reconstruction  of  the  final 
enhanced  image.  It  is  worth  mentioning  that  there  is  a  certain  overlap  between  the 
enhancement  modules:  for  example,  the  enhancement  strategy  developed  for  stellate  lesions 
may  also  enhance  circumscribed  masses  and  vice  versa  since  the  two  share  many  common 
properties. 

In  addition  to  its  efficiency,  the  algorithm  is  also  well  suited  for  further  refinements; 
optimizations  can  be  performed  for  each  type  of  malignancy  alone,  and  separately  for  the 
fusion  module. 

The  work  in  the  final  year  of  our  investigation  will  concentrate  on  extensive  testing  and 
possible  refinements  of  the  developed  algorithms. 
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